classdef shapeFunction < hgsetget
    %Clase que representa las funciones de forma a utilizar de manera
    %simbolica
    %Hereda del handle hgsetget para poder utilizar los métodos set y get
        
    properties
        
    end
    
    methods
        
        function [obj] = shapeFunction(obj)            
        end
        
        function [pg,wpg] = getGaussPointsAndWeights(obj,ndm)
            global data;
            if nargin==1
                ndm = data.mesh.ndm;
            end
            npg = data.param.npg;
            if ndm==2
                type_ele_2d = data.mesh.type;
            end
                        
            if ndm==1
                if npg==1
                    psipg = [0];
                    wpg = [2];
                elseif npg==2
                    psipg = [-1/sqrt(3) 1/sqrt(3)];
                    wpg = [1 1];
                elseif npg==3
                    psipg = [-0.7745966 0 0.7745966];
                    wpg = [5/9 8/9 5/9];
                end 
                pg = [psipg'];
            end
            if ndm==2
                if strcmp(type_ele_2d,'triangles')
                    %elemento triangular
                    if npg==1 % npg = 1
                        psipg = [1/3];
                        etapg = [1/3];
                        pg = [psipg' etapg'];
                        wpg = [1];
                    elseif npg==3 % npg = 3
                        psipg = [1/2 1/2 0];
                        etapg = [0 1/2 1/2];
                        wpg   = [1/3 1/3 1/3];
                    else
                        %por defecto npg=3
                        psipg = [1/2 1/2 0];
                        etapg = [0 1/2 1/2];
                        wpg   = [1/3 1/3 1/3];
                    end
                    pg = [psipg' etapg'];
                    wpg = wpg*1/2; %area del elemento
                end
                if strcmp(type_ele_2d,'rectangles')
                    %elemento rectangular
                    if npg==1
                        psipg = [0];
                        etapg = [0];
                        wpg = [2];
                    elseif npg==2,
                        xpg  = sqrt(3)/3;
                        psipg = [ -xpg xpg ];
                        etapg = [ 0 0 ];
                        wpg   = [ 2 2 ];
                    elseif npg==3,
                        psipg = sqrt(0.6)*[-1,0,1];
                        etapg = sqrt(0.6)*[0,0,0];
                        wpg   = (1/81)*[90,144,90];
                    elseif npg==4,
                        xpg  = sqrt(3)/3;
                        psipg = [ -xpg xpg xpg -xpg ];
                        etapg = [ -xpg -xpg xpg xpg ];
                        wpg   = [    1   1   1    1 ];
                    elseif npg==9
                        psipg = sqrt(0.6)*[-1,1,1,-1,0,1,0,-1,0];
                        etapg = sqrt(0.6)*[-1,-1,1,1,-1,0,1,0,0];
                        wpg   = (1/81)*[25*ones(1,4),40*ones(1,4),64];
                    else
                        %por defecto npg=1
                        psipg = [0];
                        etapg = [0];
                        wpg = [2];
                    end
                    pg = [psipg' etapg'];
                end
            end
        end
        
        function [dN,N,J0,detJ] = getShapeFunctions(obj,type_shfunc,xy,pgauss)
            %recibe una lista de puntos (que definen el elemento) y una
            %serie de parametros:
            %ndm: dimensiones del problema
            %type_shfunc: tipo de shape function-> 1: constante, 2:lineal, 3:cuadratica, 4: custom
            %master: 1:shfunc en coordenadas master, 0: en naturales
            %xy: lista de puntos que conforman cada elemento 
            %pgauss: puntos de gauss en donde evaluar
            % retorna las funciones de forma (N), sus derivadas (dN), el
            % jacobiano de la transformacion (J) y su determinante (detJ)
            [nen,ndm] = size(xy);
            master = 1;
            
            if ndm==1
                x1 = xy(1); 
                x2 = xy(2); 
                psi = pgauss;
                    
                if ~master
                    
                    x = (x2-x1)*(psi+1)/2; %voy a coordenadas naturales
                    
                    switch type_shfunc
                        case 1
                            %ELEMENTO CONSTANTE NATURAL
                            N(1) = 1;
                            dN(1) = 0;
                            J0 = 1;
                            detJ = 0;
                        case 2
                            %ELEMENTO LINEAL NATURAL
                            N(1) = (x2-x)/(x2-x1);
                            N(2) = (x-x1)/(x2-x1);
                            dN(1) = -1/(x2-x1);
                            dN(2) = 1/(x2-x1);
                            J0 = 1;
                            detJ = 1;
                        case 3 
                            %ELEMENTO CUADRATICO NATURAL
                            x3 = xy(:,3); 
                            J0 = 1;
                            detJ = 1;
                            
                        otherwise
                            fprintf('Caso de funcion de forma no implementado')
                    end
                else
                    switch type_shfunc
                        case 1
                            %ELEMENTO CONSTANTE MASTER
                            N(1)    = 1;
                            dN(1)   = 0;
                            J0      = 1;
                            detJ    = 0;
                            
                        case 2
                            %ELEMENTO LINEAL MASTER
                            N(1)    = (1-psi)/2;
                            N(2)    = (psi+1)/2;
                            dN1dx   = -1/2;
                            dN2dx   = 1/2;
                            dxdpsi  = x1.*dN1dx + x2.*dN2dx;
                            J0      = 1/dxdpsi;
                            detJ    = 1/dxdpsi;
                            dN = [dN1dx; dN2dx];
                        case 3
                            %ELEMENTO CUADRATICO MASTER
                             x3     = xy(3);
                             N(1)   = -0.5*(psi-psi^2);
                             N(2)   = -psi^2;
                             N(3)   = 0.5*(psi+psi^2);
                             dN1d1  = -0.5*(1-2*psi);
                             dN2d1  = -2*psi;
                             dN3d1  = +0.5*(1+2*psi);
                             dxdpsi = x1*dN1d1+x2*dN2d1+x3*dN3d1;
                             J0   = 1./dxdpsi;
                             detJ   = dxdpsi;
                             dN = [dN1d1; dN2d1; dN3d1];
                        otherwise
                            fprintf('Caso de funcion de forma no implementado')
                    end
                    
                    
                end
                
            end
            if ndm==2
                if master
                    psi = pgauss(:,1);
                    eta = pgauss(:,2);
                    
                    if type_ele_2d==1
                        switch type_shfunc
                            case 2
                                %ELEMENTO TRIANGULAR LINEAL
                                x1 = xy(1,1);
                                x2 = xy(2,1);
                                x3 = xy(3,1);
                                y1 = xy(1,2);
                                y2 = xy(2,2);
                                y3 = xy(3,2);

                                N(1) = 1 - eta - psi;
                                N(2) = eta;
                                N(3) = psi;

                                dN1d1 = -1;
                                dN1d2 = -1;

                                dN2d1 = +1;
                                dN2d2 =  0;

                                dN3d1 =  0;
                                dN3d2 = +1;   

                                dxdpsi = x1.*dN1d1 + x2.*dN2d1 + x3.*dN3d1;
                                dxdeta = x1.*dN1d2 + x2.*dN2d2 + x3.*dN3d2;

                                dydpsi = y1.*dN1d1 + y2.*dN2d1 + y3.*dN3d1;
                                dydeta = y1.*dN1d2 + y2.*dN2d2 + y3.*dN3d2;

                                J0 = [dydeta -dydpsi;
                                    -dxdeta dxdpsi];

                                detJ = dxdpsi.*dydeta-dxdeta.*dydpsi;


                                dN(1,1) = dN1d1;
                                dN(1,2) = dN1d2;
                                dN(2,1) = dN2d1;
                                dN(2,2) = dN2d2;
                                dN(3,1) = dN3d1;
                                dN(3,2) = dN3d2;

                            case 3 
                                %ELEMENTO TRIANGULAR CUADRATICO

                                x1 = xy(1,1);
                                x2 = xy(2,1);
                                x3 = xy(3,1);
                                x4 = xy(4,1);
                                x5 = xy(5,1);
                                x6 = xy(6,1);
                                y1 = xy(1,2);
                                y2 = xy(2,2);
                                y3 = xy(3,2);
                                y4 = xy(4,2);
                                y5 = xy(5,2);
                                y6 = xy(6,2);

                                N(1) = psi + eta + 2*psi^2 + 2*eta^2 + 4*eta*psi;
                                N(2) = 2*psi^2 - psi;
                                N(3) = 2*eta^2 - psi;
                                N(4) = 4*psi + 4*eta - 4*psi^2 - 4*eta*psi;
                                N(5) = 4*eta*psi;
                                N(6) = 4*psi + 4*eta - 4*eta^2 - 4*eta*psi;
                                
                                dN1d1 = 1+4*psi+4*eta;
                                dN1d2 = 1+4*psi+4*eta;
                                dN2d1 = 4*psi-1;
                                dN2d2 = 0;
                                dN3d1 = 0;
                                dN3d2 = 4*eta-1;
                                dN4d1 = 4-8*psi-4*eta;
                                dN4d2 = 4-4*psi;
                                dN5d1 = 4*eta;
                                dN5d2 = 4*psi;
                                dN6d1 = 4-4*eta;
                                dN6d2 = 4-4*psi-8*eta;


                                % Calculo del Jacobiano J0 y su determinante
                                dxdpsi = x1*dN1d1 + x2*dN2d1 + x3*dN3d1 + x4*dN4d1 + x5*dN5d1 + x6*dN6d1;
                                dxdeta = x1*dN1d2 + x2*dN2d2 + x3*dN3d2 + x4*dN4d2 + x5*dN5d2 + x6*dN6d2;
                                dydpsi = y1*dN1d1 + y2*dN2d1 + y3*dN3d1 + y4*dN4d1 + y5*dN5d1 + y6*dN6d1;
                                dydeta = y1*dN1d2 + y2*dN2d2 + y3*dN3d2 + y4*dN4d2 + y5*dN5d2 + y6*dN6d2;

                                J0 = [dydeta -dydpsi;
                                    -dxdeta dxdpsi];

                                detJ = dxdpsi.*dydeta-dxdeta.*dydpsi;

                                dN(1,1) = dN1d1;
                                dN(1,2) = dN1d2;
                                dN(2,1) = dN2d1;
                                dN(2,2) = dN2d2;
                                dN(3,1) = dN3d1;
                                dN(3,2) = dN3d2;
                                dN(4,1) = dN4d1;
                                dN(4,2) = dN4d2;
                                dN(5,1) = dN5d1;
                                dN(5,2) = dN5d2;
                                dN(6,1) = dN6d1;
                                dN(6,2) = dN6d2;
                                
                            otherwise
                                fprintf('Caso de funcion de forma no implementado')
                        end
                    else
                        switch type_shfunc
                            case 2
                                %ELEMENTO BILINEAL (RECTANGULO)
                                x1 = xy(1,1);
                                x2 = xy(2,1);
                                x3 = xy(3,1);
                                x4 = xy(4,1);
                                y1 = xy(1,2);
                                y2 = xy(2,2);
                                y3 = xy(3,2);
                                y4 = xy(4,2);
                                
                                N(1) = 1/4*(-psi-eta+eta*psi);
                                N(2) = 1/4*(psi-eta-eta*psi);
                                N(3) = 1/4*(psi+eta+eta*psi);
                                N(4) = 1/4*(-psi+eta-eta*psi);
                                
                                dN1d1 = 1/4*(-1+eta);
                                dN1d2 = 1/4*(-1+psi);
                                dN2d1 = 1/4*(1-eta);
                                dN2d2 = 1/4*(-1-psi);
                                dN3d1 = 1/4*(1+eta);
                                dN3d2 = 1/4*(1+psi);
                                dN4d1 = 1/4*(-1-eta);
                                dN4d2 = 1/4*(1-psi);


                                % Calculo del Jacobiano J0 y su determinante
                                dxdpsi = x1*dN1d1 + x2*dN2d1 + x3*dN3d1 + x4*dN4d1;
                                dxdeta = x1*dN1d2 + x2*dN2d2 + x3*dN3d2 + x4*dN4d2;
                                dydpsi = y1*dN1d1 + y2*dN2d1 + y3*dN3d1 + y4*dN4d1;
                                dydeta = y1*dN1d2 + y2*dN2d2 + y3*dN3d2 + y4*dN4d2;

                                J0 = [dydeta -dydpsi;
                                    -dxdeta dxdpsi];

                                detJ = dxdpsi.*dydeta-dxdeta.*dydpsi;

                                dN(1,1) = dN1d1;
                                dN(1,2) = dN1d2;
                                dN(2,1) = dN2d1;
                                dN(2,2) = dN2d2;
                                dN(3,1) = dN3d1;
                                dN(3,2) = dN3d2;
                                dN(4,1) = dN4d1;
                                dN(4,2) = dN4d2;
                                
                            otherwise
                                fprintf('Caso de funcion de forma no implementado')
                        
                        end
                    end                    
                else
                    fprintf('No hay implementado elementos bidimensionales en coordenadas naturales')
                    
                end
            
            end
    end
    
    end
end